function handles = excitationProfileAnalysisB(handles)
% excitationProfileAnalysisB(maps, spacingXY, somaXYZ, figFlag)
% excitationProfileAnalysisB     Analyze responses in terms of distance from soma
% Syntax:
%   excitationProfileAnalysis(maps, spacingXY, somaXYZ, figFlag)
% - maps: a map stack.
% - spacingXY is a 2-element vecor
% - somaXYZ is a 3-element vectors.
% - NB -- somaXYZ(3) is the focal plane relative to the soma plane
% - set figFlag to 1 for graphical display of map and gaussian fit
%
% Directly derived from responseDistance.m, with addition
% of stuff from analyzeArray
%
% See also responseDistance, analyzeArray

% note -- need to figure out how to analyze correctly for Z > 0
% correct way is probably
% - do the gaussian fit without factoring in Z
% - factor the Z in for the mean-distance calculation

% more notes -- 
% - could re-write distance arrays using meshgrid
% - need to revise distance arrays to accommodate non-zero soma positions

% gs Sept 2004
% gs may 2005 -- modified for mapAverager
% ********************************************************************

% get list of names of cells/expts
listOfCells = get(handles.cellList, 'String');

% type of map to analyze
% handles = getMapStr(handles);
mapStr = 'cellA.apNumber;';
handles.analysis.mapAverager.mapType = mapStr;
mapStr = handles.analysis.mapAverager.mapType;

% main loading loop
numCells = size(listOfCells, 1);
for n = 1:numCells
    
    % load the cells' data
    cellA = eval([listOfCells{n}]);     % load the cell/expt data
    mapSpacing(n) = cellA.mapSpacing;   % map spacing
    cellType{n} = cellA.cellType;       % cell type
    XY(:,:,n) = cellA.xyNew;               % soma XY position %% GS20060601 - changed to xyNew from xy
    try
        barrelXY(:,:,n) = cellA.barrelXY;   % barrel XY position
    catch
    end

%     cellA.tracesWithSpikes 
%     cellA.numberOfSpikes 
%     cellA.totalNumberOfSites 
%     cellA.totalNumberOfSpikes 
%     cellA.spikesPerSite 
%     meanDistanceFromSoma(n) = cellA.meanDistanceFromSoma; % GS20060601
%     cellA.spikeLatencies 
%     cellA.apBinary 
%     cellA.apNumber 
%     cellA.apDelay
    
    eval(['map = ' mapStr]);
%     map = cellA.genMean;                % cell's genMean map
    M(:,:,n) = map;
    
    % re-orient if necessary
    if strcmp(cellA.orientation, 'upside down') % Y values need to change sign
        %         XY(:,2,n) = -XY(:,2,n);
        %         barrelXY(:,2,n) = -barrelXY(:,2,n);
        M(:,:,n) = flipud(map);
    end
end

disp(' ');
disp('% EXCITATION PROFILE ANALYSIS:');
disp(' ');

maps = M;
spacingXY = unique(mapSpacing);
if numel(spacingXY) > 1
    disp('Error in excitation profile analysis: some cells in the group appear to have different map spacing');
    return
else
    spacingXY = [spacingXY spacingXY]; % assume for now that X and Y spacing are equal
end
somaXYZ = squeeze(XY);
if numel(unique(somaXYZ(1,:))) > 1 | numel(unique(somaXYZ(2,:))) > 1
    somaXYZ
    numel(unique(somaXYZ(1,:)))
    numel(unique(somaXYZ(2,:)))
    disp('Error in excitation profile analysis: some cells in the group appear to have different soma locations');
    return
end
somaXYZ = somaXYZ(:,1)';
if numel(somaXYZ) ==2
    somaXYZ = [somaXYZ 0];
end

figFlag = 1;

[R,C,P] = size(maps);

warning off MATLAB:divideByZero;

% DISTANCES ----------------------------------------------------------

% arrays and vectors of distances
colVals = [-((R-1)/2)*spacingXY(1):spacingXY(1):((R-1)/2)*spacingXY(1)]';
colVals = colVals - somaXYZ(1);
colArray = repmat(colVals, 1, C);

rowVals = [-((C-1)/2)*spacingXY(2):spacingXY(2):((C-1)/2)*spacingXY(2)];
rowVals = rowVals - somaXYZ(2);
rowArray = repmat(rowVals, R, 1);

planeArray = somaXYZ(3) * ones(R,C);

distArray = sqrt(rowArray.^2 + colArray.^2 + planeArray.^2);
% dists = round(distArray), return
[r,c] = size(distArray);
distances = reshape(distArray, r*c, []);
distances2 = [-distances; distances];

% RESPONSES ----------------------------------------------------------

for p = 1:P
    map = maps(:,:,p);
    [r,c] = size(map);
    responses = reshape(map, r*c, []);

    % data format for later plotting and gaussian fit
    responses1(:,p) = responses;
    responses2(:,p) = [responses; responses];

    % distances of responses > 0
    I = find(responses > 0); 
    
    % display raw data
    disp('#   R');
    disp('------');
    disp(num2str(round([responses(I) distances(I)])));
    disp('------');
    
    % mean distance weighted by response number
    resps = responses(I);
    dists = distances(I);
    distVec = [];
    for k = 1:length(I)
        distVec = [distVec; repmat(dists(k), resps(k), 1)];
    end
    
    meanWeightedDistance(p) = mean(distVec);
    stdWeightedDistance(p) = std(distVec);
    numel(distVec);
    semWeightedDistance(p) = std(distVec) / sqrt(numel(distVec));
    
    % unweighted version
    meanDistance(p) = mean(distances(I));
    stdDistance(p) = std(distances(I));
    semDistance(p) = std(distances(I)) / sqrt(numel(I));
end


% response distance analysis in terms of each unique distance

% need to re-do this -- if one map, error bars are within-map SEM; if >1 map, bars are across-cell SEM

U = unique(distances)'; % list of unique distance values
for n = 1:numel(U)     % for each unique distance ...
    I = find(distances==U(n));   % get the indices for occurrences of that distance
    % distances(I) % uncomment to verify correct indices
    A = [];
    A = responses1(I, :);
    A = reshape(A, [], 1);
%     Aavg(n) = sum(A) / P; 
    Aavg(n) = mean(A);
    Asem(n) = std(A)/sqrt(numel(A));
end
Xuniq = U;



% average map
mapAvg = mapAverage(maps, 1);

for p = 1:P
    totalNumberOfSites(p) = length(find(maps(:,:,p)));
    totalNumberOfSpikes(p) = sum(sum(maps(:,:,p)));
    normNumber(p) = totalNumberOfSpikes(p) * (spacingXY(1)) * (spacingXY(2)); % in um!!
    warning off MATLAB:divideByZero;
    spikesPerSite(p) = totalNumberOfSpikes(p) / totalNumberOfSites(p);
end

% set up figure
if figFlag
    hfig = figure(109);
    set(hfig, 'Position', [65   695   560   420], 'DoubleBuffer', 'on', 'Color', 'w');
    pltrows = 2;
    pltcols = 2;
    pltnum = 0;
    colormap jet2;
    
    % plot the average map
    pltnum  = pltnum + 1;
    hplt = subplot(pltrows, pltcols, pltnum);
    himg = imagesc(flipud(mapAvg));
    set(himg, 'XData', [min(colVals) max(colVals)], 'YData', [min(rowVals) max(rowVals)]);
    addPixInterpOption(himg);
    set(gca, 'YDir', 'normal');
    xlim([min(colVals) max(colVals)]);
    ylim([min(rowVals) max(rowVals)]);
    pbaspect([1 1 1]);
    colorbar;
    hold on;
    hsoma = plot([0 0], [0 0], 'wo');
    ylabel('Distance (um)');
    title('R-click for options');
    
    % plot distances vs responses
    pltnum  = pltnum + 1;
    hplt = subplot(pltrows, pltcols, pltnum);
    cla;
    plot(distances, responses1, 'bo');
    hold on;
    
    herr = errorbar(Xuniq, Aavg, Asem, 'k-');
    set(herr(2), 'Marker', 's', 'MarkerFaceColor', 'none', 'LineWidth', 2);
%     set(herr, 'Marker', 's', 'MarkerFaceColor', 'none', 'LineWidth', 2);
%     % for matlab7
    xlabel('Distance (um)');
    ylabel('Number of spikes');
    title('Single cell resolution');
    
    % plot folded data, gaussian fit
    pltnum  = pltnum + 1;
    hplt = subplot(pltrows, pltcols, pltnum);
    
    [r1,c1] = size(responses2);
    responses2 = reshape(responses2, r1*c1, []);
    [r2,c2] = size(distances2);
    distances2 = repmat(distances2, c1/c2, 1);

    out = responseDistanceGaussFit(hplt, distances2, responses2);
    fitting_parameters = out % for command-window display
    
    % how can you get the results in a usable data format?
    % i.e., how do you unpack the cfit class?
    % answer:
    coeffs = coeffValues(out);
    
    hfit = findobj(get(gca, 'Children'), 'Color', 'r');
    fitVals = get(hfit, 'YData');
    fitX = get(hfit, 'XData');
    peakVal = max(fitVals);
    minVal = min(fitVals);
    halfMax = (peakVal - minVal)/2;
    I = find(fitVals>halfMax);
%     figure, plot(fitVals(I), 'ro')
    fullWidth = round( fitX(I(end)) - fitX(I(1)) );
    disp('Resolution estimates from gaussian fit:');
    disp(['Full width at half-max (um): ' num2str(fullWidth)]);
    disp(['Half-width at half-max (um): ' num2str(fullWidth/2)]);

end

% analysis results
disp(' ');
disp('Resolution estimates based on mean weighted distance (um):');
disp(['Full-width at half-max (um): ' num2str(round(2*nanmean(meanWeightedDistance))) ...
        ' +/- ' num2str(round(2*nanstd(meanWeightedDistance)/sqrt(numel(meanWeightedDistance)))) ' (mean +/- sem)']);
disp(['Half-width at half-max (um): ' num2str(round(nanmean(meanWeightedDistance))) ...
        ' +/- ' num2str(round(nanstd(meanWeightedDistance)/sqrt(numel(meanWeightedDistance)))) ' (mean +/- sem)']);
disp(' ');
disp('Intensity (''brightness'') estimates based on mean weighted distance (um):');
disp(' ');

% more analysis
% disp('More analysis: ');
disp(['Spike-generating sites per map: ' num2str(nanmean(totalNumberOfSites)) ';']);
disp(['Spikes per map: ' num2str(nanmean(totalNumberOfSpikes)) ';']);
disp(['Spikes per site per map: ' num2str(nanmean(spikesPerSite)) ';']);
disp(['Spikes per map, normalized by spacing (units: AP um^2): ' num2str(nanmean(normNumber)) ';']);
excitationArea = pi * nanmean(meanWeightedDistance)^2;
disp(['Spikes per cell per flash: ' num2str(nanmean(normNumber)/excitationArea) ';']);
disp('  calculated as: (total spikes per map) * (Xspacing * Yspacing) / (pi * meanWeightedDistance^2)');
disp(' ');
disp('Note: Data stored under handles.analysis.mapAverager.ep');
disp(' ');

% output to handles
handles.analysis.mapAverager.ep.meanWtdDist_HWHM = round(10*meanWeightedDistance)/10;
handles.analysis.mapAverager.ep.totalNumberOfSites = totalNumberOfSites;
handles.analysis.mapAverager.ep.totalNumberOfSpikes = totalNumberOfSpikes;
handles.analysis.mapAverager.ep.spikesPerSite = round(100*spikesPerSite)/100;
handles.analysis.mapAverager.ep.normNumber = normNumber;
